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I study the buckling transition under compression of a two-dimensional, hexagonal, regular elastic 
honeycomb. Under isotropic compression, the system buckles to a configuration consisting of a unit 
cell containing four of the original hexagons. This buckling pattern preserves the sixfold rotational 
symmetry of the original lattice but is chiral, and can be described as a combination of three different 
elemental distortions in directions rotated 27r/3 from each other. Non-isotropic compression may 
induce patterns consisting in a single elemental distortion or a superposition of two of them. The 
numerical results compare very well with the outcome of a Landau theory of second order phase 
transitions. 



I. INTRODUCTION 

A two dimensional honeycomb structure formed by 
solid walls is the prototype of a cellular solid [1]. These 
are materials widely used in applications due to their re- 
markable mechanical properties, for instance its capacity 
for energy absorption under impact, and its low weight. 
Energy absorption is related to plastic deformation un- 
der stress. But still the ideally elastic and perfectly uni- 
form two dimensional honeycomb presents some not com- 
pletely solved puzzles. Under compressive stress, it has 
a buckling transition in which some (or all) of their walls 
bend. This transition is reminiscent of the well known 
buckling transition of an elastic bar under compressive 
stress at its extremes [2]. There has been some contro- 
versy on what the buckling mode of a regular honeycomb 
should be. On one hand, in their original work [3] , Gibson 
and Ashby presented the results of an experiment using 
an elastomeric honeycomb, under what they called biax- 
ial loading, in which they observed a non-trivial buckling 
pattern consistent with a symmetry breaking in which 
four original cells form the new repetitive motif of the 
material. In a posterior paper [4], Hutzler and Weaire 
performed numerical simulations and did not observe this 
pattern, but instead a buckling mode equivalent to that 
obtained under uniaxial stress. They argue that the pat- 
tern observed in Ref. [3] was a consequence of finite 
size effect, and the use of flat confining walls. Numer- 
ical results taking into account these effects [5] did show 
the pattern observed by Gibson and Ashby. Very re- 
cently, Okomura, Ohno and Noguchi [6] have studied the 
problem using a combination of a homogenization tech- 
nique and finite elements numerical simulations. Their 
results do not agree with those of Hutzler and Weaire 
[4]. Instead, they found buckling patterns that can be 
interpreted as a superposition of three individual buck- 
ling modes. They also found that whether one, two, or 
three of these modes are active depends on the degree of 
anisotropy of the externally applied strain. 

In view of the aforementioned contradiction between 



[4] and [6] , and considering that the techniques employed 
in both cases are quite different, an independent inves- 
tigation to determine which of the two results is correct 
seems appropriate. In the first (numerical) part of this 
paper I will show that appropriately done numerical sim- 
ulations using the technique used in [4] do not support 
the results claimed there, but instead those reported in 
[6] . In the second (more theoretical) part I will show how 
the results obtained in the simulations arc fully compat- 
ible with the predictions of a Laundau theory of second 
order phase transition applied to the buckling problem. 
This theory allows to obtain at once the buckled config- 
uration of the system under a generic form of the macro- 
scopically homogeneous applied deformation. 

II. NUMERICAL SIMULATION 

I have simulated a two dimensional honeycomb 
through the technique used in Refs. [4], [5], namely by 
considering the honeycomb walls as one dimensional rods, 
and including stretching and bending energy as 

1 f ( oil \ 
E stretch = -k s J dl 

E b end = l;h J c 2 dl (1) 

where c is the local curvature. To discretize these ex- 
pressions I have used 7 intermediate points between any 
two neighbor vertices, but particular cases where checked 
using 18 points, to guarantee the absence of noticeable 
effects due to discretization. The only essential param- 
eter of the model is kt/L 2 k s , where L is the length of 
the individual rods. This ratio is physically related to 
the fraction A of two dimensional space that is occupied 
by the rods [1,4], namely A = 4^/fcf,/ 'L 2 k s . The simula- 
tions presented below were done at kb/L 2 k s — 4.5 x 10 -4 
(then A = 0.085), but additional checks indicate that 
the results are not qualitatively dependent on the precise 
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value of the parameter. All previously obtained buckling 
patterns for perfect honeycombs can be accommodated 
within a 2 x 4 unit cell. Then the elemental cell I simu- 
lated is precisely the 2x4 cell shown in Fig. 1(a), with 
periodic boundary conditions. The simulation method 
consists in calculating the forces acting on all points of 
the discretized system, and updating their positions us- 
ing a viscous dynamics. The control variable was the 
macroscopic strain, that can be changed varying the size 
and shape of the simulation box. Stresses in the sys- 
tem can be evaluated both by numerical differentiation 
of the total energy with respect to strain, and by direct 
summation in terms of the forces between particles. The 
equivalence of the two results allows to check for consis- 
tency and convergence of the simulation. 

Before indicating the results obtained, it is clarifying 
to discuss qualitatively the behavior observed (see [6]). 
The buckling structures that appear are related to reac- 
commodation of the vertices of the honeycomb structure, 
in such a way that lines of vertices forming zig-zag chains 
shift relatively to neighbor chains, as qualitatively indi- 
cated in Fig. 1(b). There are three of these modes, that 
will be refereed to as the elementary modes of buckling. 
The patterns they generate will be called the uniaxial 
patterns. They are characterized by the unitary vectors 
shown in Fig. 1(e). Whether one, two, or the three ele- 
mentary modes acquire non-zero amplitude at the buck- 
ling transition depends on the macroscopic strain ap- 
plied. The qualitative pictures of Fig. 1(c) and (d) show 
the effect of combining two or three elementary modes. 
We see that the patterns of Fig. 1(b) and (c) are iden- 
tical (except for the wall bending, that in this qualita- 
tive picture is not taken into account) to the patterns in 
Ref. [3]. The pattern in Fig. 1(c) will be referred to 
as the Gibson- Ashby pattern. The pattern in Fig. 1(d), 
that combines the three elementary modes, preserves the 
hexagonal symmetry of the structure (this will be referred 
to as the symmetric pattern). A buckling mode of this 
symmetry has been observed and simulated in [7], but 
only in the case of plastic buckling. To my knowledge, 
the only prediction of this mode for a perfectly elastic 
honeycomb is contained in [6] . It is remarkable that this 
pattern has lost the mirror symmetry plane of the orig- 
inal structure: it is a chiral pattern. The chirality can 
be defined as the sign of the product of the amplitudes 
of the three elementary modes from which the pattern is 
constructed. 

In the numerical simulations, starting from the un- 
strained structure and isotropically compressing it, I have 
observed that the symmetric pattern is the stable con- 
figuration after buckling. However, it has to be men- 
tioned that the transition from the unbuckled to the 
symmetrically buckled configuration, although continu- 
ous (then implying no surmounting of any energy barrier) 
may take a long computational time, and that other pat- 
terns (which correspond to saddles of the energy) can be 



transitory observed. For instance, I show in Fig. 2 the 
mechanical energy of the system as a function of the com- 
pressive strain s in a single run increasing s at a constant 
rate. At the buckling point (at s ~ 0.0037, note that I 
define s as s = (s x + s y )/2), and in the particular run 
shown, the system seems to buckle to the Gibson- Ashby 
pattern shown in Fig 3(b). This is not a minimum but 
a saddle of the energy, however it lasts long enough, and 
the energy of this branch can be numerically followed (as 
seen in Fig. 2 with open triangles) up to some strain at 
which it is observed to transform to the real energy min- 
imum which corresponds to the symmetric pattern. In 
fact, when preparing the system in the symmetric pat- 
tern for large strain, and upon reduction of the strain, we 
obtain the results indicated by full circles in Fig. 2, which 
correspond to the real energy minimum of the system. In 
addition, if the system is prepared in the uniaxial or in 
the Gibson-Ashby pattern at large strain, and strain is 
reduced (sufficiently rapidly for the configuration not to 
destabilize, and sufficiently slowly so as energy can be 
calculated accurately) we can follow the energy of these 
two configurations, as indicated by the stars, and the full 
triangles in Fig. 2. 

The uniaxial and Gibson-Ashby patterns may be- 
come the stable buckling modes under appropriate non- 
isotropic loading. For instance, by compressing along the 
x (y) direction, and keeping the perpendicular direction 
unstrained, I have observed the Gibson-Ashby (uniaxial) 
pattern to be stable after buckling. The stability of the 
Gibson Ashby pattern under particular uniaxial loading 
agrees again with the results in [6], but not with those 
claimed in [4] . We will see now how all these results can 
be fully systematized within a theoretical framework. 

III. LANDAU THEORY OF THE BUCKLING 
TRANSITION 

The results presented in the previous section are 
enough numerical input to construct the Landau theory 
of this peculiar second order transition. In fact, we see 
in Fig. 2 that the energies of all buckled configurations 
(whether the true minimum or the saddles) meet together 
in value and derivative at the buckling point, coinciding 
also at that point in value and derivative with the branch 
corresponding to the unbuckled system. This means that 
at the buckling point, in a generic parameter space, the 
state point of the system passes from a configuration with 
a single minimum (for the unbuckled state) to one with 
different minima and saddles in a continuous manner. 

I will present a Landau description in which the order 
parameters are the three (small) amplitudes fa, fa, fa 
of the elementary modes. For convenience, these three 
modes will be associated with the unitary vectors v 1 , v 2 , 
and v 3 shown in Fig. 1(e). The free energy of the system 
(in the present case it actually corresponds simply to the 
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elastic energy) must a scalar, and then it can only con- 
tain combinations of the amplitudes and the (eventually 
anisotropic) external strains that are invariant with re- 
spect to the symmetry operations of the lattice. Consid- 
ering for the moment only the isotropically compressed 
case, we should look for invariant combinations of the 
amplitudes. Up to fourth order those available are: 
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Then in this case the most general form of the free energy 
describing a second order transition is 

F = a(s c - s) (0 2 + <& + 4) + (3 (0? + 4 + cftf + 

-7(^! + 0i0l + ^?) (3) 



with numerical constants a, (3, and 7. This free energy 
has a single minimum at (j) p = (p = 1,2,3) for s < 
s c , representing the unbuckled state. For s > s c , this 
expression has saddles (or minima) at the following values 
of the amplitudes: 



= A, 
= B, 



= 0, 
= B, 



b r = {) 
6 r = 



4>l = C, cj) q = C, <f> r = c 



(4) 
(5) 
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where p, q, r are arbitrary permutations of 1, 2, 3, and 
A, B, C, are given by: 
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(s - s c )a 
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B = 


(s - s c )a 
4/3-7 


(8) 


C = 


(s - s c )a 
6/3-27 


(9) 



They correspond respectively to the uniaxial, Gibson- 
Ashby, and symmetric patterns. The corresponding val- 
ues of the free energy are, 
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(I have subtracted the free energy of the unbuckled sys- 
tem f un b, that is taken as zero in the Landau theory, but 
that should be included when comparing with the results 
of the numerical simulations). We see that the symmetric 
pattern is the minimum energy one for 7 > (whereas 
the uniaxial pattern provides the absolute minimum if 



7 < 0). Then, in order to agree with the simulation re- 
sults, we will assume 7 > 0. The physical justification for 
the positive sign of 7 is not provided by the Landau the- 
ory, and should come from an explicit evaluation of the 
mechanical energy of the honeycomb. From the previous 
values of the free energy a parameter free relation can be 
obtained and compared with the numerical results. First 
note that he ratio (3/j can be obtained for instance as 
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This in particular should be independent of the actual 
value of the compression s (as long as s > s c ). From the 
numerical results it can be obtained that this ratio is in 
fact independent on s, and the numerical value is approx- 
imately 7//3 ~ 0.01 for the parameters of the simulation. 
Then note that to a very good approximation we have 
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This is a parameter free relation that has to be satis- 
fied in our numerical simulations. From the data in Fig. 
2, it can be in fact verified that this is very accurately 
satisfied. This is a strong evidence that the present Lan- 
dau theory describes the physics of the present buckling 
transition. 

In order to make the theory more complete, I want 
to consider now the possibility of non-isotropic external 
loading on the system. This means that instead of a 
single parameter s, we have now a generic (symmetric) 
strain tensor (i, j — 1, 2) applied onto the system (the 
previously introduced isotropic compression s is related 
to the trace of this tensor). This has to be introduced 
into the free energy in a symmetrically invariant form. 
To lowest order I will include it only in the second order 
term, which is the one that triggers the transition. Using 
the unitary vectors v 1 , v 2 , vu 3 of the elementary modes, 
two different terms quadratic in the amplitudes can be 
written: 



i,j=l,2p,g=l,2,3 
i,j=l,2p,9=l,2,3 



pq'rp'rq 



pqVp'+'q 



(15) 
(16) 



where a pq and b pq are arbitrary numeric matrices. How- 
ever, these expressions have to be invariant under permu- 
tation of the elementary vectors (since this is a symme- 
try operation obtained by the mirror symmetry along the 
line containing the third vector) and sign change of any of 
the amplitudes (which is obtained by a particular spatial 
translation allowed by symmetry). Then it is obtained 
that both a pq and b vq matrices should be proportional to 
the identity, i.c, 
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(17) 
(18) 



Then becomes proportional to Su{4>\ + (fy + </>§) 
and is the term considered in the isotropically compressed 
case. 

To analyze the second contribution it can be more con- 
venient to use the following definition of the three inde- 
pendent components of the strain tensor: 

s = (sn + s 22 )/2 

52 = (sil - «22)/2 

53 = S12 = S21 

which represent the applied deformation in a more phys- 
ical way: s represents an isotropic compression (we al- 
ready used this), whereas s 2 and S3 are the two indepen- 
dent shear modes, which are related by a 7r/4 rotation. 
In terms of these variables, and using explicitly the com- 
ponents of the unitary vectors we finally arrive to the 
following form of the free energy: 



a [(s c - s) (<\>\ + 4>l + <f>t) 
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s 2 {<& - 4>j/2 - 4/2) + s 3 ^- (4 - <g) 



+ 



+ p {<t>\ + 4>l + tl) 2 - ^\<g + fyl + <f>l<fi) (19) 

This is the final expression for the free energy close to 
the buckling transition. Minimizing it we can obtain the 
buckled state under any particular combination of the 
three independent strains s, s 2 , and S3. Note that as this 
expression has to describe also the unstrained system, 
and then as c and S will be respectively proportional to 
bulk and shear modulus of the original structure, then 
typically S << as c . 

I want to describe now the buckling mode map of the 
system, namely, what the amplitudes of the three ele- 
mentary modes are for any choice of the strain tensor. 
First notice the following scaling of the free energy: If 
we consider the values of the three order parameters at 
the minimum of (19), namely <^ lin , to be a function of 
s c — s, s 2 , and S3, then the following relation is satisfied: 

Cm( S c - S, S 2 , S3) = \~ 1/2 (f> P mm (HSc - S), As 2 , As 3 ) 

(20) 

This implies in particular that the borders between dif- 
ferent regions in the parameters space s c — s, s 2 , and S3 
are spanned by rays propagating from the origin. 

I will analyze a couple of particular cases. First con- 
sider the case s 3 = 0, i.e., purely compressive strains 
along x and y (although non necessarily equal). I show 
in Fig. 4(a) the map of buckling modes in the s-s 2 plane 



for this case. The borders between different regions can 
be worked out analytically. All of them are straight lines 
emanating from the point s = s c , s 2 = 0, as the previous 
argument indicates. The transition between unbuckled 
{4>. p = 0) and uniaxial pattern {<j>\ ^ 0) can be easily ob- 
tained setting (f>2 — </>3 — in (19). The limit line is given 
by s 2 = (s — s c )a/S. The transition line between the un- 
buckled and the Gibson- Ashby pattern (0 2 = 7^ 0) is 
obtained along the same lines, as s 2 = — 2(s — s c )a/S. 
Increasing s at s 2 = 0, the symmetric pattern appears 
at s = s c , as we already know from the isotropically 
compressed case. The symmetric pattern looses its strict 
rotational symmetry for any s 2 7^ 0. However, it still 
has the three elementary model active as long as we are 
within the V-shaped region in Fig. 4(a). The limits of 
this region are given by 
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for the transition to the uniaxial pattern, and 



1 cry s c — s 
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(21) 



(22) 



for the transition to the Gibson- Ashby pattern. Note the 
exact relation 



(GA — >sym) / (uni- 
A 7*2 



= -2 



(23) 



valid for any values of the parameters of the free energy. 

The results of Fig. 4(a) are fully compatible with the 
results in [6] (in particular, relation (23) is very well sat- 
isfied). We note that for the present parameters the sta- 
bility of the single elementary mode and Gibson-Ashby 
pattern that was obtain numerically under uniaxial com- 
pression is recovered. 

As an additional example I show the map of buckling 
modes in the s 2 -S3 plane for some s > s c in Fig. 4(b). 
Note the nice symmetry of this pattern, which has one 
two or three elementary modes active depending on the 
particular choice of the applied strains s 2 and S3 (re- 
member that s 2 and S3 are related by a rotation of 7r/4). 
Again, all borders between different sectors are straight 
lines. The analytical expression for the line separating 
sectors 1 and 1,2 is given by 



S3 



cry 
73/3<5 



(S c - S) 



(24) 



All other lines can be obtained then from symmetry. 

To finish, we note that all transitions in the parameter 
space are continuous, namely, there are no jumps of the 
order parameters at any point, and there is no possibility 
of metastabilities either. 
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IV. CONCLUSIONS 



The buckling mode of an elastic two dimensional hon- 
eycomb provides an example of non-trivial patterns with 
symmetry breaking appearing in a very simple mechan- 
ical system. Remarkably, for isotropic compression the 
symmetry breaking produces the appearance of a chiral 
ground state. This problem is also a realization of a 
second order transition that can be accurately modeled 
through a Landau theory constructed on the basis of the 
symmetry of the problem. The agreement between the 
Landau theory and the numerical simulation is seen to 
be very good. 
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FIG. 1. (a) The hexagonal starting lattice. Dotted box is 
the system actually simulated, (b) Upon shifts of the vertices 
as indicated by the arrows, the uniaxial pattern is obtained. 
Note however that there are three equivalent ways of gener- 
ating this pattern, that can be characterized by the unitary 
vectors shown in (e). Combining two or the three of them we 
obtain the configurations in (c) and (d). 
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FIG. 2. Total energy of the simulated system (dotted 
box in Fig. 1(a)) as a function of compressive strain 
s = (s x + s y )/2. A global, linear contribution in s has been 
subtracted to appreciate tiny differences in energy. All con- 
tinuous lines are quadratic fitting of the date close to the 
buckling point. Results are shown from a running from the 
unbuckled state (squares), in which a Gibson- Ashby pattern 
is generated first (open triangles). This is however a saddle 
of the energy, and is seen to transform to the symmetric pat- 
tern after some time. The other three curves show the results 
starting from the uniaxial (stars), Gibson- Ashby (full trian- 
gles), and symmetric (circles) patterns at high strain, and 
reducing it towards the unstrained configuration. Note how 
the minimum energy of the buckled state is obtained for the 
symmetric pattern. The letters indicate where the snapshots 
in Fig. 3 were taken. 
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FIG. 4. The buckling modes map in the S-S2 plane for 
S3 = (a), and in the S2-S3 plane for a constant value of s > s c 
(b). In each region, the numbers indicate which elementary 
modes are active (see Fig. 1). The analytical expressions for 
the limits between different regions are given in the text. 



FIG. 3. The uniaxial, Gibson- Ashby, and symmetric buck- 
ling modes. The snapshots correspond to the points indicated 
in Fig. 2. For isotropic compression, the symmetric pattern 
provides the minimum energy, the other two correspond to 
saddles, and eventually destabilize. However, they can be 
made stable under appropriate non-isotropic loading (the dis- 
placements with respect to the hexagonal configuration have 
been amplified by a factor of 5 to render the geometrical struc- 
ture more visible). 
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